Getting started

Adapted from: SOSTOOLS' SOSDEMO1 (See Section 4.1 of SOSTOOLS User's Manual) and Example 2.4 of [PJ08]

P. Parrilo and A. Jadbabaie Approximation of the joint spectral radius using sum of squares. Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402

using DynamicPolynomials
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

\[ 5y^{4} - x^{2}y^{2} + 2x^{3}y + 2x^{4} \]

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
con_ref = @constraint(model, p >= 0)
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We see above that the solver found a feasible solution. We now inspect this solution:

q = gram_matrix(con_ref)

GramMatrix with row/column basis:
 MonomialBasis([y^2, x*y, x^2])
And entries in a 3×3 SymMatrix{Float64}:
  5.0                     4.4408920985006264e-17  -2.9518518518518517
  4.4408920985006264e-17  4.903703703703703        1.0
 -2.9518518518518517      1.0                      2.0

We can get the SOS decomposition from the gram matrix as follows:

sosdec = SOSDecomposition(q)

(-2.1236350925148186*y^2 + 0.685602733872361*x*y + 1.4043287476933908*x^2)^2 + (-0.6931517747987985*y^2 - 2.1053486943394355*x*y - 0.020343251432301056*x^2)^2 + (0.09856272587979267*y^2 - 0.03405099490006257*x*y + 0.16567112157245334*x^2)^2

We now seek for the SOS decomposition of the following polynomial:

p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

\[ y^{2} + x^{2} - xy^{2} + 4x^{4}y^{6} \]

We build the same model as previously with this new polynomial. Here we can use Model instead of SOSModel as we explicitly constrain p to belong to the SOS cone with p in SOSCone().

model = Model(solver)
con_ref = @constraint(model, p in SOSCone())
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can query the SOS decomposition directly from the constraint reference as follows:

sos_decomposition(con_ref)

(-0.8633544111907666*y + 0.1929775871076749*x*y + 2.811891400614671e-18*x*y^2 + 1.974804063246114*x^2*y^3)^2 + (-1.0044005280851495e-16*y + 0.7962855107812749*x - 3.604670114761173e-17*x*y - 1.7452642425911866*x*y^2 - 5.2795306922881863e-17*x^2*y^3)^2 + (0.24383463418983614*y - 1.7639164207923265e-16*x - 1.548639212897575*x*y + 3.447822046557963e-17*x*y^2 + 0.2579336224358902*x^2*y^3)^2 + (-5.040390840687567e-17*y + 0.604920974441955*x + 4.554029379037467e-17*x*y + 0.27599821010522485*x*y^2 - 4.6874560662018034e-17*x^2*y^3)^2 + (-0.44177350740730575*y - 5.425371220674706e-17*x - 0.10009637589813167*x*y + 1.3534821192250128e-17*x*y^2 - 0.18335527863613868*x^2*y^3)^2

This page was generated using Literate.jl.